
/** @file polyhedra.cc
 *
 * Implementation of polyhedra 
 *
 */

/*
 *  Copyright (C) 2007 Stefan Weinzierl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include "polyhedra.h"
#include "constants.h"

namespace sector_decomposition {

  // -----------------------------------------------------------------------------------

  /**
   *
   * Constructor.
   *
   */
  polyhedra::polyhedra() :
    index(), coordinates()
  {}

  /**
   *
   * Constructor.
   *
   * Initialises the polyhedron from the polynomial poly with the help of arg_index and global_data.
   *
   */
  polyhedra::polyhedra(GiNaC::ex poly, const std::vector<size_t> & arg_index, const integration_data & global_data) :
    index(arg_index.size()), coordinates()
  {
    size_t n = arg_index.size();

    // initialize with 0,1,...,n-1
    for (int k=0; k<n; k++) index[k]=k;

    poly = poly.expand();

    // only one monomial ?
    if ( !GiNaC::is_a<GiNaC::add>(poly) ) 
      {
	std::vector<GiNaC::ex> point;

	GiNaC::ex monom = poly;

	for (size_t k=0; k<n; k++)
	  {
	    int deg = monom.degree(global_data.list_feynman_parameter[arg_index[k]]);

	    point.push_back(GiNaC::numeric(deg));
	  } // k=0; k<n; k++

	coordinates.push_back(point);
      }
    // now we can assume that poly has more than one term
    else 
      {
	size_t j_poly = poly.nops();

	for (int j=0; j<j_poly; j++)
	  {
	    std::vector<GiNaC::ex> point;

	    GiNaC::ex monom = poly.op(j);

	    for (size_t k=0; k<n; k++)
	      {
		int deg = monom.degree(global_data.list_feynman_parameter[arg_index[k]]);

		point.push_back(GiNaC::numeric(deg));
	      } // k=0; k<n; k++

	    coordinates.push_back(point);
	  } // j
      } // more than one term
  }

  /**
   *
   * Merges the two polyhedra p1 and p2 into one (new) polyhedra.
   *
   */
  polyhedra::polyhedra(const polyhedra & p1, const polyhedra & p2) :
    index(p1.index), coordinates()
  {
    for (int j=0; j<p1.coordinates.size(); j++)
      {
	coordinates.push_back(p1.coordinates[j]);
      }

    for (int j=0; j<p2.coordinates.size(); j++)
      {
	coordinates.push_back(p2.coordinates[j]);
      }
  }

  /**
   *
   * Returns the number of points stored in the array coordinates.
   */
  size_t polyhedra::number_corners() const
  {
    return coordinates.size();
  }

  /**
   *
   * Removes redundant points from the coordinate list.
   *
   * Redundant points are points which are not corners.
   *
   */
  polyhedra & polyhedra::keep_corners_only()
  {
    size_t n = index.size();

    size_t n_points = coordinates.size();

    std::vector<bool> flag_remove(n_points, bool(false));

    for (int j1=0; j1<n_points; j1++)
      {
	const std::vector<GiNaC::ex> & point_1 = coordinates[j1];

	for (int j2=j1+1; j2<n_points; j2++)
	  {
	    const std::vector<GiNaC::ex> & point_2 = coordinates[j2];

	    if ( (!flag_remove[j1]) && (!flag_remove[j2])  )
	      {
		GiNaC::ex local_d_max = 0;
		GiNaC::ex local_d_min = 0;

		for (size_t k=0; k<n; k++)
		  {
		    GiNaC::ex d1 = point_1[k];
		    GiNaC::ex d2 = point_2[k];

		    if (d1-d2 > local_d_max)
		      {
			local_d_max = d1-d2;
		      }
		    else if (d1-d2 < local_d_min)
		      {
			local_d_min = d1-d2;
		      }
		  } // k=0; k<n; k++
		
		if ( local_d_min >= 0 )
		  {
		    flag_remove[j1] = true;
		  }
		else if ( local_d_max <= 0 )
		  {
		    flag_remove[j2] = true;
		  }

	      } // !flag_remove
	  } // j2
      } // j1

    for (int j=n_points-1; j>=0; j--)
      {
	if ( flag_remove[j] ) coordinates.erase( coordinates.begin()+j );
      }

    return *this;
  }

  /**
   *
   * Shifts the polyhedra, such that in the corresponding polynomial
   * no single variable can be factorized.
   *
   * In Spivakovsky's notation this corresponds to
   * \f[
   *    \Delta \rightarrow \tilde{\Delta}.
   * \f]
   *
   */
  polyhedra & polyhedra::factorize_single_variables()
  {
    size_t n = index.size();

    size_t n_points = coordinates.size();

    // find omega
    std::vector<GiNaC::ex> omega = coordinates[0];

    for (int j=1; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = coordinates[j];

	for (size_t k=0; k<n; k++)
	  {
	    if ( point[k] < omega[k] ) omega[k] = point[k];
	  } // k
      } // j

    // shift polyhedra
    for (int j=0; j<n_points; j++)
      {
	std::vector<GiNaC::ex> & point = coordinates[j];

	for (size_t k=0; k<n; k++) point[k] -= omega[k];
      } // j

    return *this;
  }

  /**
   *
   * Adds a point to the polyhedron.
   *
   */
  polyhedra & polyhedra::add_point(const std::vector<GiNaC::ex> & point)
  {
    coordinates.push_back(point);

    return *this;
  }

  /**
   *
   * Returns the factorization degree.
   *
   */
  GiNaC::ex polyhedra::get_factorization_degree() const
  {
    GiNaC::ex f_degree = 0;
    bool flag_ini = true;

    size_t n = index.size();

    size_t n_points = coordinates.size();

    if ( n_points == 0 ) return Infinity;

    for (int j=0; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = coordinates[j];

	GiNaC::ex local_degree = 0;
	for (size_t k=0; k<n; k++) local_degree += point[k];

	if ( (local_degree<f_degree) || flag_ini )
	  {
	    f_degree = local_degree;
	    flag_ini = false;
	  }
      } // j

    return f_degree;
  }

  /**
   *
   * Returns the factorization degree for the set index_set.
   *
   */
  GiNaC::ex polyhedra::get_factorization_degree(const multi_index_ordered & index_set) const
  {
    GiNaC::ex f_degree = 0;
    bool flag_ini = true;

    size_t n = index_set.size();

    size_t n_points = coordinates.size();

    for (int j=0; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = coordinates[j];

	GiNaC::ex local_degree = 0;
	for (size_t k=0; k<n; k++) local_degree += point[index_set[k]];

	if ( (local_degree<f_degree) || flag_ini )
	  {
	    f_degree = local_degree;
	    flag_ini = false;
	  }
      } // j

    return f_degree;
  }

  /**
   *
   * Determines the vector omega.
   *
   */
  void polyhedra::get_omega(std::vector<GiNaC::ex> & omega) const
  {
    size_t n = index.size();

    size_t n_points = coordinates.size();

    // find omega
    omega = coordinates[0];

    for (int j=1; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = coordinates[j];

	for (size_t k=0; k<n; k++)
	  {
	    if ( point[k] < omega[k] ) omega[k] = point[k];
	  } // k
      } // j
  }

  /**
   *
   * Fills a vector S of boolean variables, whether or not a variable 
   * occurs in monomials of degree f_degree.
   *
   */
  void polyhedra::find_S(std::vector<bool> & S, GiNaC::ex f_degree) const
  {
    size_t n = index.size();

    size_t n_points = coordinates.size();

    S.clear();
    for (size_t k=0; k<n; k++) S.push_back(false);

    for (int j=0; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = coordinates[j];

	GiNaC::ex local_degree = 0;
	for (size_t k=0; k<n; k++) local_degree += point[k];

	if ( local_degree==f_degree )
	  {
	    for (size_t k=0; k<n; k++) 
	      {
		if ( point[k] > 0 ) S[k]=true;
	      }
	  }
      } // j
  }

  /**
   *
   * The projection.
   *
   */
  polyhedra & polyhedra::projection(const std::vector<bool> & S, GiNaC::ex f, GiNaC::ex c)
  {
    size_t n = index.size();

    size_t n_points = coordinates.size();

    // store old variables
    std::vector<size_t> temp_index = index;

    std::vector<std::vector<GiNaC::ex> > temp_coordinates = coordinates;

    index.clear();
    coordinates.clear();

    // set index
    for (size_t k=0; k<n; k++)
      {
	if ( !(S[k]) ) index.push_back( temp_index[k] ); 
      }

    // set coordinates
    for (int j=0; j<n_points; j++)
      {
	const std::vector<GiNaC::ex> & point = temp_coordinates[j];

	GiNaC::ex local_S_degree = 0;
	for (size_t k=0; k<n; k++) 
	  {
	    if ( S[k] ) local_S_degree += point[k];
	  }

	if ( local_S_degree < f )
	  {
	    std::vector<GiNaC::ex> new_point;

	    for (size_t k=0; k<n; k++) 
	      {
		if ( !(S[k]) ) new_point.push_back( point[k]*c/(f-local_S_degree) );
	      }

	    coordinates.push_back(new_point);
	  } // local_S_degree < f

      } // j

    return *this;
  }

  /**
   *
   * Finds a minimal subsector, such that index_set is permissible.
   *
   */
  void polyhedra::find_minimal_subsector(std::vector<size_t> & index_set, GiNaC::ex f_degree) const
  {
    size_t n = index.size();

    bool flag_continue = true;

    for (int m=1; m<=n; m++)
      {
	multi_index_ordered local_multi_index(n,m);

	for( local_multi_index.init(); !local_multi_index.overflow(); local_multi_index++)
	  {
	    if ( flag_continue )
	      {
		if ( get_factorization_degree(local_multi_index) >= f_degree )
		  {
		    for (int k=0; k<m; k++) index_set.push_back( local_multi_index[k] );

		    flag_continue = false;
		  }
	      }
	  } // loop over local_multi_index
      } // m=1; m<=n;
  }

  /**
   *
   * Output. Prints out the index vector and the coordinates.
   *
   */
  std::ostream & operator<< (std::ostream & os, const polyhedra & arg)
  {
    os << "class polyhedra" << std::endl;

    os << " Index: ";
    for (int j=0; j<arg.index.size(); j++) os << arg.index[j] << " ";
    os << std::endl;

    os << " Coordinates: " << std::endl;
    for (int j=0; j<arg.coordinates.size(); j++) 
      {
	const std::vector<GiNaC::ex> & point = arg.coordinates[j];
	for (int k=0; k<point.size(); k++) os << " " << point[k];
	os << std::endl;
      }

    return os;
  }

} // namespace sector_decomposition

